Fractionation effects in phase equilibria of polydisperse hard sphere colloids 
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The equilibrium phase behaviour of hard spheres with size polydispersity is studied theoretically. 
We solve numerically the exact phase equilibrium equations that result from accurate free energy 
expressions for the fluid and solid phases, while accounting fully for size fractionation between co- 
existing phases. Fluids up to the largest polydispersities that we can study (around 14%) can phase 
separate by splitting off a solid with a much narrower size distribution. This shows that exper- 
imentally observed terminal polydispersities above which phase separation no longer occurs must 
be due to non-equilibrium effects. We find no evidence of re-entrant melting; instead, sufficiently 
compressed solids phase separate into two or more solid phases. Under appropriate conditions, 
coexistence of multiple solids with a fluid phase is also predicted. The solids have smaller poly- 
dispersities than the parent phase as expected, while the reverse is true for the fluid phase, which 
contains predominantly smaller particles but also residual amounts of the larger ones. The proper- 
ties of the coexisting phases are studied in detail; mean diameter, polydispersity and volume fraction 
of the phases all reveal marked fractionation. We also propose a method for constructing quantities 
that optimally distinguish between the coexisting phases, using Principal Component Analysis in 
the space of density distributions. We conclude by comparing our predictions to perturbative theo- 
ries for near-monodisperse systems and to Monte Carlo simulations at imposed chemical potential 
distribution, and find excellent agreement. 

PACS numbers: 82.70.Dd, 64.10.+h, 82.70.-y, 05.20.-y 
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I. INTRODUCTION 
A. The hard sphere model 

Hard spheres are particles that do not interact except 
via an infinite repulsion on contact. In a hard sphere 
system there is no contribution to the internal energy, 
U , from interparticle forces since U is zero for all the 
allowed configurations. Minimising the free energy, F = 
U — TS, is thus equivalent to maximising the entropy, 
S: the structure and phase behaviour of hard spheres 
is determined solely by entropy. Temperature T only 
features as a trivial factor setting the energy scale. 

The hard sphere model was originally introduced as 
a mathematically simple model of atomic liquids (see 
e.g. but has since also been recognised as a useful ba- 
sic model for complex fluids Q such as spherical colloids. 
Colloidal particles coated with a thin polymeric layer so 
that strong steric repulsions dominate the attractive dis- 
persion forces between the colloidal cores behave in many 
ways as hard spheres. Indeed, crystallisation can be ob- 
served at densities similar to those predicted by computer 
simulation for hard spheres, with a single-phase fluid be- 
low volume fractions <f> « 0.494, fluid-solid coexistence at 
up to <j> 0.545, and a single-phase solid at higher volume 
fractions Measurements of the osmotic pressure 

and compressibility similarly show very good agreement 
with predicted hard sphere properties p|. 

There is, however, one important and unavoidable dif- 



ference between colloids and the classical hard sphere 
model: whereas the spheres in the classical model are 
identically sized, colloidal particles have an inevitable 
spread of diameters. The magnitude of this spread is con- 
veniently characterised by the parameter 8, which is often 
also referred to as polydispersity and measures the stan- 
dard deviation of the diameter distribution normalised 
by its mean: 



(1) 



Here the averages a and a 2 are defined via 



pa 1 = J dap(a)a t 



(2) 
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with pier) the density distribution of the system. The lat- 
ter is defined so that the number density of particles with 
diameter between a and a + da is given by p(a) da. The 
total density is then p = J dap(a), and n(a) — p(a)/p 
is the normalised diameter distribution. The /?,; are the 
moments of the density distribution, with po = p. The 
presence of polydispersity in the system brings in a new 
parameter that allows us to distinguish between size dis- 
tributions of different widths; the shape of the diameter 
distribution is of course also relevant. Compared to the 
monodisperse case, polydispersity causes several qualita- 
tively new phenomena which have received much interest 
in recent years. 
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FIG. 1: Sketch of the fluid (F) and solid (S) phase boundaries 
for polydisperse hard spheres, following |18j . The boundaries 
are plotted as polydispersity 8 versus volume fraction <j>. The 
fluid boundary approaches the solid one until they meet at a 
terminal polydispersity, S t . For S just below S t , this scenario 
suggests re-entrant melting: compressing the crystal to suffi- 
ciently high volume fraction should transform it back into a 
fluid. 



B. New phenomena arising from polydispersity 

The effect of polydispersity on the phase behaviour of 
hard spheres has been investigated by experiments 0, 0] , 
computer simulations 0, M- 13. lift ITl| , density functional 
theorie s Il2l flU and simplified analytical theories [ToL 
Ellllllilllllllllli3- We will now outline the main 
findings and introduce the relevant terminology. 

First, it is intuitively clear 01 that significant diam- 
eter polydispersity should destabilize the crystal phase, 
because it is difficult to accommodate a range of diam- 
eters in a lattice structure. Experiments have indeed 
shown that crystallisation is suppressed above a termi- 
nal polydispersity of 8t ~ 0.12 [3, @. Since then much 
theoretical work has focused on estimating <5 t . Dickinson 
et al 0, for example, extrapolated the decrease of the 
volume change on melting with polydispersity to zero, 
obtaining an estimate of St ~ 0.12. Pusey 0] used 
a simple Lindemann-type criterion to estimate that the 
larger spheres in a polydisperse system would disrupt the 
crystal structure above S t ~ 0.06... 0.12. McRae and 
Haymet used density functional theory (DFT) and 
found that there was no crystallisation above 5 t ~ 0.05. 
Barrat and Hansen ^2] also employed DFT, estimating 
the free energy difference between fluid and solid. Taken 
together, this body of theoretical work suggests that the 
terminal polydispersity arises from a progressive narrow- 
ing of the fluid-solid coexistence region with increasing 6, 
with phase boundaries meeting at St |l 'il Il5| in a point 
that has been identified as one of equal concentration 0] 
(rather than a critical point). 

Bartlctt and Warren 18] also found re-entrant melting 
on the high-density side of this point: for 8 just below 
St , they predicted that compressing a crystal could trans- 
form it back into a fluid. Fig. ^ shows a sketch of this 
scenario. 

Physically, the existence of re-entrant melting would 
suggest that, while in the monodisperse case the solid 



has the lower free energy at all volume fractions above 
4> ~ 55%, the fluid can become preferred again at large 
</> if the polydispersity is sufficiently large. This result 
is compatible with the intuition that polydispersity re- 
duces the maximum packing fraction in a crystal (since 
a range of diameters need to be accommodated on uni- 
formly spaced lattice sites), while it increases the maxi- 
mum packing fraction in the fluid, where smaller spheres 
should be able to fill "holes" between larger particles 
more easily. 

This intuition can be made more quantitative by com- 
paring the fluid and solid free energies, following [2l| . 
The basic analysis by Bartlett and Warren 0] ignores 
fractionation, i.e. the fact that coexisting phases need 
not have identical diameter distributions as long as they 
combine to give the correct overall or "parent" distri- 
bution //°)(er). The normalised diameter distribution is 
thus fixed and equal in all phases. In the moment free 
energy (MFE) method described below this corresponds 
to retaining only the overall density p . Phase equilibria 
can then be found by the usual double-tangent construc- 
tion |20j from a plot of the (moment) free energy density 
/ versus po- We display such free energy plots in Fig. El 
showing along the cc-axis the volume fraction <p rather 
than pq; the two are proportional for fixed diameter dis- 
tribution. (The free energies are those also used for our 
detailed calculations below. The diameter distribution 
was of a Schultz form, but other size distributions are 
expected to give similar results.) For polydispersities up 
to S = 0.08 a tangent plane between the solid and fluid 
phases always exists, giving the conventional fluid-solid 
coexistence. At S = 0.08, re-entrant melting has ap- 
peared: a second double tangent is possible because at 
large volume fractions the solid free energy is now higher 
than that of the fluid. As S increases, the solid free energy 
continues to increase relative to the fluid and eventually 
lies above the latter for all (see S = 0.09 and 0.1 in 
Fig. [2J • The point where this first happens gives the ter- 
minal polydispersity 8 t ; in our example S t ~ 0.083. As 
5 approaches St from below, the widths of both the or- 
dinary and the re-entrant fluid-solid coexistence regions 
shrink to zero and merge into the point of equal con- 
centration, consistent with the phase diagram shown in 
Fig.Q] 

As mentioned above, this picture ignores the possi- 
bility of fractionation. Bartlett and Warren inves- 
tigated fractionation effects approximately, by using a 
MFE with two density variables included, po and p\ — 
pa. They concluded that the phase diagram topology 
remained qualitatively unchanged; quantitatively, the 
point of equal concentration was shifted to higher den- 
sity and lower polydispersity. It has to be born in mind, 
however, that while the approach of [0| allowed coexist- 
ing phases to have different mean diameters, it implicitly 
still constrained them to have the same 8. (This is be- 
cause, within the MFE method applied to the Schultz 
prior R(a) cx <7 z e~( z+1 ' <T of [la, the density distribu- 
tions p(o-) = R(a)e Xo+Xia oc o-*e\ Xl ~( z+1)]a in all phases 
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FIG. 2: Free energy / versus 4> within the "no fractionation" 
approximation used by Bartlett and Warren [hH| . Phase sep- 
aration occurs where double tangents between the fluid (thick 
line) and solid (dashed line) branches of the free energy can 
be drawn; the plot at 8 — 0.04 shows an example (thin line). 
At 8 = 0.08 re-entrant melting can be observed: two double 
tangents can now be drawn. For larger 8, phase separation is 
no longer predicted. (A linear term —190 has been added to 
all free energies to make the plots more readable.) 



are also Schultz, with common z and therefore common 
S = (1 + z)^ 1 / 2 .) On the other hand, numerical simula- 
tions that allow for fractionation show that a solid with 
a narrow size distribution can coexist with an essentially 
arbitrarily polydisperse fluid H, 0, |nj. This suggests 
that the prediction of re-entrant melting should be re- 
examined theoretically, allowing for such fractionation 
effects. Conceptually, it also implies that the concept 
of a terminal polydispersity is likely to be useful only 
for the solid but not for the fluid, and we will see this 
confirmed below. 

Fractionation has also been predicted to lead to solid- 
solid coexistence pH Il7l l2l| . where a broad diameter 
distribution is split into a number of narrower solid frac- 
tions. This occurs because the loss of entropy of mix- 
ing is outweighed by the better packing, and therefore 
higher entropy, of crystals with narrow size distribution; 
accordingly, as the overall polydispersity of the system 
grows, the number of coexisting solids is predicted to in- 
crease. Fig.|21sketches this effect, following the treatment 
of ^(| . There is no coexistence region between fluid and 
solid, due to a simplification in the analysis of [l6l | : rather 
than solving the phase equilibrium conditions, only the 
free energies were equated between the fluid and the (one 
or several) solid phases. The resulting lines in the phase 
diagram generally lie inside the actual phase separation 
region, but give a rough guide to the phase transitions 
that can occur. The parent diameter distribution consid- 
ered had a "top hat" form (uniform between given min- 
imum and maximum diameters), and for multiple solids 
fractionation was assumed to be "hard" , with the parent 
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FIG. 3: Sketch of fluid (F) and multiple solid (S) phase coex- 
istences in polydisperse hard spheres, following jig. Approx- 
imate phase boundaries are plotted as polydispersity 8 versus 
volume fraction <j>. For sufficiently large 8 and <j>, coexistence 
of several solids is predicted; see text for discussion. 



distribution split into non-overlapping top hat distribu- 
tions with identical polydispersities. 

Previous work as described above leaves open a number 
of questions. The rather drastic, and differing, approxi- 
mations for size fractionation used in previous studies of 
re-entrant melting and solid-solid coexistence [lit UtI ITsL 
l2lj . as described above, leave the relative importance of 
these two phenomena unclear. Theoretical calculations 
that account fully for fractionation remain restricted to 
highly simplified van der Waals free energies [ijj. Nu- 
merical simulations can in principle also capture arbi- 
trary fractionation behaviour, but have been carried out 
at constant chemical potential distribution 0, 0, ^| . As 
explained in more detail in Section IVI Bl the system's 
overall particle size distribution can then change dramat- 
ically across the phase diagram. This is in contrast to the 
experimental situation and so limits the applicability of 
the results. 

Our main aim in this study is, therefore, to calcu- 
late the equilibrium phase behaviour of polydisperse hard 
spheres on the basis of accurate free energy expressions, 
taking full account of fractionation and going beyond 
previous work on fluid-solid and solid-solid coexistence. 
The experimentally observed behaviour of hard sphere 
colloids will of course also depend on non- equilibrium ef- 
fects, e.g. the presence of a kinetic glass transition |22| . 
anomalously large nucleation barriers |2^] or the growth 
kinetics of polydisperse crystals (24|. Nevertheless, the 
equilibrium phase behaviour needs to be understood as a 
baseline from which non-equilibrium effects can be prop- 
erly attributed. Also, more of the equilibrium behaviour 
may be observable under microgravity conditions, where 
the glass transition is shifted to higher densities or even 
absent [2q . 

We begin in Section [H] by defining the free energies we 
use to describe the fluid and solid phases of polydisperse 
hard spheres. Section llHl reviews the moment free energy 
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TABLE I: Relations between dimensional and dimensionless 
quantities. All dimensional quantities except for the units (3, 
o"o and wo themselves are denoted by tildes "~" . 



method and its numerical implementation for solving the 
phase equilibrium conditions. In Section llVI we then de- 
scribe the basic features of the phase behaviour that we 
find; a short account of these results has appeared in [2^| . 
Section[V]describes in detail the fractionation effects that 
we predict, and introduces a new method for constructing 
optimal visualisations of polydisperse phase behaviour. 
Section Ivl! finally, compares our results to perturbative 
theories for the near-monodisperse limit and to Monte 
Carlo simulations at constant chemical potential differ- 
ences. The agreement is very good, thus validating our 
approach. We conclude in Section IVIII with a summary 
and outlook towards future work. 



II. FREE ENERGIES 

Our starting point is the decomposition of the free en- 
ergy of a polydisperse system into an ideal and an excess 
part, 



dap(a) [lnp(a)-l]+r*({ Pi }) 



(3) 



The excess part / cx can, in principle, depend on all de- 
tails of p(a) and therefore on all of its moments pi, but 
we will be concerned with truncatable free energies |27j . 
For these, the dependence is only through a finite number 
of moments, for us specifically pQ,...,p3. 

Strictly speaking, equation 10 gives the free energy 
density; we will continue to refer to this as the free energy 
for short. Also, all quantities in are dimensionless: 
we assume that sphere diameters are measured in units 
of some reference value a®, that all densities are made 
dimensionless by multiplying by the volume vq — 7tctq/6 
of a reference sphere, and that all energies are measured 
in units of T = Boltzman's constant fee is set to 1 
throughout. Free energy and pressure are then in units 
of T/vo, for example. Table [I] summarises the relations 
between important dimensionless and dimensional quan- 
tities. Conveniently, with our choice of units p% = 4> is 
simply the volume fraction of spheres. 

For the fluid phase of polydisperse hard spheres, the 
most accurate free energy approximation available is the 
generalisation by Salacuse and Stell |28( of the equation of 



state due to Boublik, Mansoori, Carnahan, Starling and 
Leland (BMCSL) ,2G| 3CJ; for the monodisperse case this 
reproduces the Carnahan-Starling equation of state [3l| . 
In our dimensionless quantities, the BMCSL expression 
for the excess free energy takes the form 



r 



4 
pi 



- p ln(l - p 3 ) 



pI 



3P1P2 
1-P3 P3(l-P3) : 



(4) 



As anticipated above, this is truncatable, involving only 
the moments pi = J dap(a)a l (i = . . . 3) of the density 
distribution. Bartlett [32j | provided an elegant argument 
why — at least within a virial expansion — such a moment 
structure of the excess free energy for the hard sphere 
fluid should in fact be exact. 

For phase coexistence calculations we will also need to 
have a compact expression for the excess free energy of 
the polydisperse hard sphere crystal. This is not at all 
a trivial question. In principle, the structure of a poly- 
disperse crystal could be rather complex, with different 
sites inside the crystalline unit cell occupied preferen- 
tially by particles with different ranges of diameters. The 
system would then effectively be an ordered solid solu- 
tion (see e.g. H3,|jM|]). Most theoretical work makes the 
simplifying assumption that one has a substitutional^ 
disordered solid, where crystal sites are assumed to be 
occupied equally likely by particles of any diameter (see 

e.g. mM)- 

A simple-minded but popular approach to estimating 
the free energy is cell theory, first introduced by Kirk- 
wood jH3| and widely used since (see e.g. 0|) : particles 
are treated as independent but confined to an effective 
cell formed by their neighbours. However, it is clear that 
for a polydisperse system this is unlikely to be a use- 
ful approximation. For example, the cells of the model 
would have to be made large enough to accommodate the 
particles with the largest diameter, even if the fraction 
of such particles is very small. 

We follow instead the more quantitative, "geometric" 
approach proposed by Bartlett |32| ■ He assumed that 
the excess free energy of the solid depends on the same 
moments po, ■ ■ ■ , P3 as that of the fluid. This can be moti- 
vated from scaled particle theory £3, > which suggests 
that the excess chemical potential, /x ex (cr), of spheres of 
diameter a is given by a cubic polynomial in a 



p 2 a 



p 3 a 



(5) 



The coefficients /ig x and /if x can be determined from 
the Widom insertion principle |^}. The latter can be 
stated as saying that exp(— pf x (a)) is the ratio of the 
(excess parts of the) partition functions for N + 1 and 
N particles, where the added particle has diameter a. 
(Equivalently the excess chemical potential may be in- 
terpreted as the work of inserting an (N + l)-th hard 
sphere of diameter a into a system of N spheres.) In a 
system with purely hard interactions, this implies that 
/i GX (cr) is positive and an increasing function of a. For 
large a, the presence of the added particle effectively 



5 



just reduces the volume available to the N others, giv- 
ing p ex (a) s» Per 3 (= P(ir/6)a 3 in dimensional units), 
hence /z| x = P. For small a, one notes that the ratio of 
the (excess) partition functions is also the average Boltz- 
mann factor of the added particle, the average being over 
the Boltzmann distribution of the TV-particle system. In 
the hard sphere case, exp(— /z ex (cr)) is thus the probabil- 
ity of being able to insert a particle without overlap. In 
the limit of vanishing particle this probability is 1 — <f), 
giving p c *(a — ► 0) = /iQ X = - ln(l — 4>). 

One now notes that |J5J implies that the excess free en- 
ergy can only depend on the moments po, ■ ■ ■ , Pz- Indeed, 
from the definition of the excess chemical potentials and 
with the dependence of the excess free energy on p(a) 
expressed through a (possibly infinite) set of moments 
Pi, 



Sf cyi 



df c 



dpi 



A comparison with the form (JSJ of the excess chemical 
potentials reveals that / ex can only depend on po, . . . , P3, 
as claimed. The same is then true also for the /i° x — 
df ex /dpi (i = 0, . . . , 3), which are recognised as excess 
moment chemical potentials. The excess free energy of 
a polydisperse hard sphere mixture can thus be deduced 
from that of any other mixture which is equivalent in 
the sense of having the same po, ■ ■ ■ , ps- These moments 
determine the number density along with the basic geo- 
metric properties of mean particle diameter, surface area 
and volume. The simplest mixture with a finite num- 
ber of species that can match any given po , ■ ■ ■ , P3 is a 
bidisperse one. Indeed, this has four degrees of freedom, 
namely, the number densities and particle diameters of 
the two species. We can therefore identify the excess free 
energy of a polydisperse hard sphere solid with that of 
the equivalent bidisperse system. For the latter, we fol- 
low Bartlett in using the fits to the simulation data of 
Kranendonk et al (36J. Because these data are obtained 
for an fee substitutionally disordered crystal, an implicit 
assumption is that the polydisperse crystal will have the 
same structure. 

There is a difficulty in Bartlett 's approach with the 
determination of the excess moment chemical potentials 
^q x , . . . , pf^. He fixed ^q x and /x| x to the exact re- 
sults derived from the Widom insertion principle, /ig x = 
— ln(l — pa) and pJ^ = P. The remaining two excess mo- 
ment chemical potentials, /xf x and /i§ x , can then be found 
from the bidisperse simulation data, by requiring /i ex (cr) 
at the diameters of the small and large spheres to agree 
with the simulated excess chemical potentials of the two 
species. However, because of the approximate character 
of the excess free energy, the p ex (cr) derived by this route 
do not obey the thermodynamic consistency requirement 
5p e *(a)/5p(a') = Sp ex (a')/Sp(a), which corresponds to 
dpi*/ dp j = dp° x /dpi for the excess moment chemical po- 
tentials. To avoid this in our study, we assign the latter 
by explicitly evaluating the derivatives of the excess free 
energy, p1 x = df c */dpi. Thermodynamic consistency is 



then automatic. The price we pay is that our p ex (a) 
no longer has the theoretically expected asymptotic be- 
haviour for (7 — > and a — * oo. This means that we have 
to restrict use of our solid free energy to relatively nar- 
row diameter distributions, as discussed in more detail 
below. 



III. NUMERICAL METHOD 

A. Moment free energy 

Our computational approach for determining the phase 
behaviour of polydisperse hard spheres is based on the 
moment free energy (MFE) method. We giv e a brief 
outline here; details can be found in [2(3, H3, El E3- 
Recall first the phase equilibrium conditions for coexis- 
tence of p phases, in a system described by a truncat- 
able free energy. By definition, the excess free energy 
then depends on a finite set of M (generalised) moments 
Pi = J da p(a)wi(a) defined by weight functions Wi(cr); 
above we had u>i{a) = a 1 . In coexisting phases, the chem- 
ical potentials p(a) and pressure P must be equal. The 
former are, by differentiation of ©, 



p{a) 



Sf 
5p(a) 



lnp(<r) +^/xfwi(cr) 



(G) 



with /i° x = <9/ ox / dpi as before. The pressure is given by 
the Gibbs-Duhem relation 

P=-f+ [da p(a)p(a) = Po - / cx + ]T p? P i (7) 



To the conditions of equality of chemical potentials and 
pressure we need to add the requirement of conservation 
of particle number for each species a, which reads 



E 



(8) 



where a = 1, . . . ,p labels the phases and is the frac- 
tion of the system volume occupied by phase a. One 
then finds from equality of the p(<j), Eq. ©, together 
with particle conservation JHJ), that the density distribu- 
tions in coexisting phases can be written as 



Here the x[ a ^ must obey 



exp 






£ 7 f (7) exp 





(9) 



A 



(a) 



-//,- 



(«).• 



+ Cj 



(10) 



and the Ci are undetermined constants that do not affect 
the density distributions (0- One can fix them e.g. by re- 
quiring all the \[ a ^ in one of the phases to be zero. A little 
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reflection then shows that l|10|) together with J2 a v = 1 
and the equality of the pressures Q) in all phases give a 
closed system of nonlinear equations for the p(M + 1) 

variables X^ and v^ a \ A solution can thus, in princi- 
ple, be found by a standard algorithm such as Newton- 
Raphson. Generating an initial point from which such 
an algorithm will converge, however, is still a nontrivial 
problem, especially when more than two phases coexist 
and/or many moments pi are involved. Furthermore, the 
nonlinear phase equilibrium equations permit no simple 
geometrical interpretation or qualitative insight akin to 
the construction of phase diagrams from the free energy 
surface of a finite mixture. 

The moment free energy addresses these two disadvan- 
tages. To construct it, one starts by modifying the free 
energy decomposition @ to 



/ = J da p(a) 



In 44-1 



+ /° x ({pJ) 



(11) 



In the first (ideal) term, a normalising factor R(a) has 
been included inside the logarithm. This has no effect on 
the exact thermodynamics because it contributes only 
terms linear in p(a), but will play a central role below. 
One can now argue that the most important moments 
to treat correctly in the calculation of phase equilibria 
are those that actually appear in the excess free energy 
/° x ({Pi})- Accordingly, one imposes particle conserva- 
tion JSJ only for the pi, but allows it to be violated in 
other details of the density distribution p(a) which do 
not affect the pi . These "transverse" degrees of freedom 
are instead chosen to minimise the free energy (jlljl . and 
more precisely its ideal part since the excess contribution 
is a constant for fixed values of the pi . This minimisation 
gives 



p(a) = R(a) exp 



^ XiWi(a) 



(12) 



where the Lagrange multipliers Xi are chosen to give the 
desired values of the moments 



Pi 



da uii (a) R(a) exp 



(13) 



The corresponding minimum value of / as given in 111|) 
then defines the moment free energy (MFE) 



/mom({ft}) 



(^2 XtPt 



po + r({ Pl }) (w) 



Since the Lagrange multipliers are (at least implicitly) 
functions of the moments pi, the MFE depends only on 
the pi. These can now be viewed as densities of "quasi- 
species" of particles, allowing for example the calculation 
of moment chemical potentials [27j 



Pi 



dpi 



Xi 



dp 



dpi 



Xi 



Pi 



(15) 



and the corresponding pressure P — J^iPiPi — /mom 
which turns out to be identical to the exact expres- 
sion (0. A finite-dimensional phase diagram can thus 
be constructed from / mom according to the usual tan- 
gency plane rules, ignoring the underlying polydisperse 
nature of the system. Obviously, though, the results now 
depend on R(a). To understand its influence, one notes 
that the MFE is simply the free energy of phases in which 
the density distributions p(a) are of the form (|12fl . To 
ensure that the parent phase is contained in the fam- 
ily, one normally chooses its density distribution as the 
prior, R(a) = p (0) (ct); the MFE procedure will then be 
exactly valid whenever the density distributions actually 
arising in the various coexisting phases are members of 
the corresponding family 



p{a) — p^(a) exp 



X t Wi{a) 



(16) 



It is easy to show from © that this condition holds 
whenever all but one of a set of coexisting phases are 
of infinitesimal volume compared to the majority phase. 
Accordingly, the MFE yields exactly the onset of phase 
of coexistence, conventionally represented via cloud and 
shadow curves (see below). Similarly, one can show that 
spinodals and critical points are found exactly |27| . 

For coexistences involving finite amounts of different 
phases the MFE only gives approximate results, since 
different density distributions from the family (TSJ, cor- 
responding to two (or more) phases arising from the same 
parent p^ (a), do not in general add to recover the parent 
distribution itself. Moreover, from Gibbs' phase rule, a 
MFE depending on M moments will not predict more 
than M + 1 coexisting phases, while we know that a 
polydisperse system can in principle separate into an ar- 
bitrary number of phases. Both of these shortcomings 
can be overcome by including extra moments within the 
MFE. By choosing the weight functions of the extra mo- 
ments adaptively, the properties of the coexisting phases 
can then be predicted with in principle arbitrary accu- 
racy H3,E3|- Importantly for us, the results can in fact 
be used as initial points from which a solution of the 
exact phase equilibrium problem can be converged suc- 
cessfully [44l |45| . This is the technique that we use here. 
Once a phase split for a given parent distribution p^ (a) 
has been found, care needs to be taken to check that it 
is globally stable, i.e. that no phase split of lower free 
energy exists [27J ■ Adopting this procedure, we are able 
to calculate coexistence of up to five phases, which so 
far has been possible only for much simpler free energies 
depending on a single density moment (see e.g. [27|L 



B. Implementation 

We focus below on parent distributions with unit mean 
particle diameter a; any other choice could be absorbed 
into the unit length a§. For small polydispersity 5, the 
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standard moments p, = Jda p{a)(J l then become very 
close to each other, and in fact strictly identical in the 
limit 8 — > 0. This causes numerical difficulties, and 
we therefore work instead with the centred moments 
p\ = J da p(a)[{a — l)/<5o] 4 which remain distinct even 
for small 8. The factor Sq is included to ensure that the 
moments are all of comparable magnitude. We therefore 
choose it in the middle of the range of polydispersities 
5 that we study, with typically So = 0.05. The centred 
moments are obviously linearly related to the conven- 
tional ones, e.g. pi = (p% — po)/So- The BMCSL and 
solid free energies can therefore readily be re-expressed 
in term of the centred moments. Because the transfor- 
mation between the two sets of moments is linear, the 
corresponding sets of excess moment chemical potentials 
Pi* = df CK /dpi are also linearly related and easily con- 
verted into each other. 

We combine the fluid and solid branches of our excess 
free energy by simply taking the minimum for a given 
set of moments. Some care is needed here: because the 
solid free energy is derived from fits to simulation data for 
bidisperse systems (see above ), we expect it to be reliable 
only in the region spanned by the simulations [36]. The 
smallest diameter ratio investigated in the simulations 
is a — 0.85. The maximum polydispersity that can be 
reached in a bidisperse system for this diameter ratio is 
5 = (a + 1/a - 2) 1 / 2 /2 ~ 0.08. We therefore restrict 
use of the solid free energy to phases with polydispersity 
below this value. Reassuringly, we will see below that 
all solid phases occurring in equilibrium phase splits are 
well below this threshold. 

A further constraint on the use of the solid free en- 
ergy arises from the fact that, as explained above, our 
excess chemical potentials p c *(cr) do not have the correct 
limiting behaviour predicted from the Widom insertion 
principle for a — > and er — > oo. Physically, this is again 
plausible because we are extrapolating to sphere diame- 
ters far from the mean of the distribution, and therefore 
far from the regime where the simulation data will be 
reliable. We will therefore always work with diameter 
distributions with hard cutoffs either side of the mean 
so that the behaviour of p ex (<r) for very small or large 
<7 never comes into play. Finally we also restrict the 
volume fractions for the solid branch of the free energy 
to 0.494 < 4> < 0.74, which are respectively the small- 
est and largest for which monodisperse hard spheres at 
equilibrium exhibit a crystalline solid phase. 



IV. PHASE BEHAVIOUR 

We now describe our results for the overall phase be- 
haviour of polydisperse hard spheres. Our numerical 
work requires a choice to be made for the parental di- 
ameter distribution. We focus mostly on a triangular 



distribution, where n^\a) = p^(a)/ p is given by 

r,(°)(Vi = — / a ~ ~ w "> for 1 ~ w - a - 1 

y > w 2 \ (1 + w) - a for l<cr<l + w 

whose width parameter w is related to the polydisper- 
sity by w — >/6 8. For the moderate values of 8 of 
interest here one expects other distribution shapes to 
give qualitatively similar results, based on the intuition 
that for narrow size distributions 5 is the key param- 
eter controlling the phase behaviour |14| . To verify 
this, we also consider briefly a Schultz parent distribu- 
tion, n(°)(cr) cx cr z e~( z+1 ) cr , cut off outside the range 
a e [0.8,1.2]. For a narrow distribution, i.e. large z, 
where the cutoffs are unimportant, the polydispersity is 
then related to the parameter z by S 2 = l/(z + 1) and 
the mean diameter is unity as before. 

A. Onset of phase coexistence 

The most basic question we can ask about phase be- 
haviour regards the onset of phase separation coming 
from single-phase regions. Increasing the volume frac- 
tion (j) of the parent at given polydispersity 8, a single- 
phase fluid (F) will first separate into coexisting fluid 
and solid (S) phases at the so-called cloud point. The lo- 
cus of all cloud points in the (0,i5)-plane defines the fluid 
cloud curve. The incipient phase corresponding to the 
cloud point is called the "shadow" solid; its properties 
define the solid shadow curve. These curves are shown 
in Fig.0](a) for a triangular parent distribution. An im- 
portant feature is that the fluid cloud curve continues 
throughout the whole range of polydispersities that we 
can investigate numerically: even at 8 = 0.14, a hard 
sphere fluid will eventually split off a solid on compres- 
sion. This is in marked contrast to the phase diagram 
of pj| as sketched in Fig. ^ The key difference is that 
our analysis accounts fully for fractionation: Fig.0|shows 
that the coexisting shadow solid always has a relatively 
modest polydispersity with 8 never rising above 0.06 
even when the cloud fluid has 8 = 0.14. This fractiona- 
tion effect prevents the convergence of the solid and fluid 
phase boundaries which a theory with fixed 8 0] pre- 
dicts, along with the resulting re-entrant melting (Fig-QJ. 
These findings are in qualitative accord with Monte Carlo 
simulations for the simpler case of fixed chemical poten- 
tials H, |jj O] , discussed in detail in Section IVIBI be- 
low. The results imply, in particular, that the terminal 
polydispersity St cannot be defined as the point beyond 
which a fluid at equilibrium will no longer phase sepa- 
rate; 8t only makes sense as the maximum polydispersity 
at which a single solid phase can exist (see below). 

The fractionation effects described above can be seen 
more explicitly by comparing the (normalised) diame- 
ter distributions of the fluid cloud and solid shadow 
phases, as displayed in Fig.[3]for (parental) polydispersity 
8 = 0.05 and 8 = 0.10. At the cloud point the size distri- 
bution in the fluid coincides with the parent distribution 
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FIG. 4: Cloud curves (thick) and shadow curves (thin), for 
polydisperse hard spheres with a triangular (a) and Schultz 
(b) diameter distribution. The curves show polydispersity 8 
versus volume fraction <j> for the cloud and shadow phases; 
dashed lines link sample cloud-shadow pairs. The solid (S) 
cloud curve has two branches, with onset of F-S and S-S coex- 
istence at low and high volume fractions, respectively. Where 
they meet, a triple point occurs; squares mark the cloud phase 
and the two coexisting shadows there. In the Schultz plot, 
the dotted lines indicate the expected continuations of the 
fluid cloud and corresponding shadow curve beyond the re- 
gion where our numerical methods work reliably. 

as it must. The distribution in the coexisting shadow 
solid, on the other hand, deviates increasingly from that 
of the parent as the parental 5 increases. In particular, 
the solid contains predominantly the larger particles and 
has a rather more narrow spread of sizes, consistent with 
the small solid polydispersities found above. We will see 
shortly that these properties are rather generic and per- 
sist inside the coexistence region. 

We now assess the effect of the shape of the particle size 
distribution on these results. Fig. 0] shows that the fluid 
cloud and solid shadow curves are qualitatively and even 
quantitatively very similar for the triangular and Schultz 
distributions. (Numerically, we can only reach 5 = 0.10 



FIG. 5: Normalised size distributions n(a) = p(er)/po for the 
coexisting fluid and shadow phases at the fluid cloud point, for 
parent distributions with polydispersities S = 0.05 (top) and 
S — 0.10 (bottom) and triangular (left) and Schultz (right) 
shape. Solid lines show the cloud fluid, which is identical to 
the parent, and dashed lines the shadow solid. 



for the latter, but have no reason to expect that this is a 
physical feature and indicate the expected continuation 
of the curves by dotted lines.) Fig. [S] (right) demon- 
strates that the qualitative features of the fractionation 
behaviour are also the same between the two distribu- 
tions, consistent with our intuition that variations in the 
shape of the parental size distribution have, for given 5, 
only a minor effect. 

We next consider the onset of phase separation coming 
from the single-phase solid, which defines the solid cloud 
curve and corresponding shadow curve. Initially we focus 
again on the triangular size distribution. Fig. 01 (a) shows 
that a decrease in density at low polydispersities leads to 
conventional fluid-solid phase separation. At higher 5, 
however, the solid cloud curve acquires a second branch 
at higher densities. This is broadly analogous to the re- 
entrant phase boundary found in 18|, but with the cru- 
cial difference that the system phase separates into two 
solids rather than a solid and a fluid. The two branches 
meet at a triple point. Here the solid cloud phase coex- 
ists with two shadow phases, one fluid and one solid, as 
marked by the squares in Fig.0|(a). The triple point is 
located at 6 ~ 0.07; since it is at the maximum of both 
branches of the solid cloud curve, this value also gives 
the terminal polydispersity beyond which solids with tri- 
angular diameter distribution are unstable against phase 
separation. 

Fig. (left) displays the diameter distributions for 
the cloud and shadow solids, at S = 0.05 on the high- 
density branch of the solid cloud curve. In comparison 
with Fig. |SJ what is striking is that the fractionation 
effects at the onset of solid-solid coexistence are much 
stronger than for fluid-solid phase separation at the same 
S. This is consistent with physical intuition. The fluid- 




FIG. 6: Normalised size distributions n(a) of solid cloud and 
shadow phases, on the high-density branch of the solid cloud 
curve at 8 = 0.05. Left: triangular parent; right: Schultz 
parent. The solid lines show the cloud phase, the dashed 
lines the shadow. Note the strong size fractionation effects. 



solid phase separation exists even in the monodisperse 
limit. The presence of polydispersity acts as a small per- 
turbation to this transition, certainly at low 5, so that 
fractionation effects can be viewed as incidental. Solid- 
solid phase separation, on the other hand, is driven by 
polydispersity and could not take place without fraction- 
ation. 

We compare again at this stage with the results for 
the Schultz parent distribution. Fig. 01 (b) shows that 
the cloud and shadow curves look qualitatively similar 
to the triangular case. Quantitatively, the low-density 
branch of the solid cloud curve now has a maximum, 
giving the terminal polydispersity as 6 t ~ 0.06. The 
triple point is at slightly smaller 5, and the whole high- 
density branch of the solid cloud curve - which describes 
the onset of solid-solid phase separation - is shifted to 
smaller S compared to the triangular parent case. Fig. 
(right) shows the diameter distributions for the cloud 
and shadow solids, at the onset of phase separation at 
6 = 0.05. Compared to the triangular parent, the frac- 
tionation effects are now even stronger. In fact, the size 
distribution of the shadow solid continues to increase to- 
wards larger sphere diameters a and is terminated only 
by the hard cutoff at a — 1.2 which we impose in the 
Schultz case; note that in the cloud solid (solid line) this 
cutoff is hardly discernible. In the triangular case, there 
is no sharp cutoff effect on the shadow solid: the form 
of the parent forces all size distributions to drop to zero 
continuously at the upper end. 

The above observations for the Schultz parent suggest 
an analogy with recent results for isotropic-nematic phase 
separation in hard rod- like particles 0, ^(| . For suffi- 
ciently wide rod length distributions, one observes there 
that the shadow nematic phase can become dominated 
by the longest rods in the system, i.e. those with lengths 
near the cutoff, even though these make up only a small 
fraction of the parent distribution. Such cutoff effects are 
important only near the cloud point: as soon as the new 
phase occupies a nonzero fraction of the overall system 
volume, particle conservation prevents it from containing 
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FIG. 7: Fractional system volume occupied by the newly 
forming solid, v^ 2 \ versus the parent volume fraction, <j>, for 
cutoffs imposing three different ranges of particle sizes a as 
indicated in the legend. Once enough of the new solid phase 
exists (above tr 2 ' ~ 0.08) the behaviour is essentially cutoff- 
independent. The cloud point, on the other hand, where v' 2 ' 
extrapolates to zero, depends strongly on the cutoff; this is 
more clearly visible in the linear- log plot in the inset. 



an atypically large number of long rods. To test whether 
we have a similar situation for the onset of solid-solid sep- 
aration from a Schultz parent, we have varied the cutoff 
on the sphere diameters. Fig. [7| plots the fractional sys- 
tem volume i/ 2 -* occupied by the new solid against the 
parent colloid volume fraction. We observe that v^ 2 > is 
indeed cutoff-independent well inside the coexistence re- 
gion, where it is non-negligible. The position of the cloud 
point itself, on the other hand, where extrapolates to 
zero, is strongly cutoff-dependent. We conclude, there- 
fore, that at the onset of solid-solid coexistence from a 
Schultz parent with S — 0.05 the shadow solid is cutoff- 
dominated, in analogy with the shadow nematics in the 
Onsager model of long hard rods 0, 0] . 

One may wonder whether the cutoff effects described 
above are an artefact of the approximate nature of our 
excess chemical potentials for the polydisperse solid. We 
cannot give a definitive answer to this question here, but 
suggest that such effects might in fact be expected. From 
(|6*|l , equilibrium of chemical potentials between the cloud 
(parent) solid p(°' (a) and the shadow solid p^ (a) implies 



p W(a) = (0 (° ) ( ( r)exp(-A/i c 



(17) 



with p ex (cr) = 2j=n A/^er 1 and A/i° x the differences 
in the excess moment chemical potentials between the 
shadow and cloud phases. The Widom insertion princi- 
ple tells us that ^| x , being equal to the pressure, is equal 
in the two phases. Thus A/i GX (u) should generically be a 
quadratic polynomial in a. If the cr 2 -term has a negative 
coefficient, then from <|17|) it will overwhelm the exponen- 
tial tail of the Schultz parent p^{a). The shadow den- 
sity distribution p^ (a) then increases strongly at large 
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a so that its properties will be dominated by the pres- 
ence of any cutoff. In fact, this argument suggests that 
the same could happen even with a (sufficiently poly- 
disperse) Gaussian parent. Only a stronger decay, say 
p(°)(<7) ^ exp(— a a ) with a > 2, could definitely prevent 
cutoff effects on the shadow solid. This question deserves 
further study, but would require more accurate excess 
chemical potentials for polydisperse solids - and over a 
larger range of sphere diameters - than we currently have 
at our disposal. 

So far we have investigated the global stability of sin- 
gle phases, i.e. the stability against macroscopic phase 
separation. One can also ask about local stability of 
the phases, i.e. the spinodal points. Since our free en- 
ergy is spliced together from separate fluid and solid 
branches, we cannot investigate instability to fluid-solid 
separation. The stability of single-phase fluids against 
fluid-fluid demixing has been studied by Warren [42| and 
Cuesta They found, using the BMCSL free energy, 
that spinodal instabilities do indeed occur, but only for 
very broad diameter distributions such as log-normals 
with 5 above « 2.5, or bimodal distributions with a wide 
size disparity between the larger and smaller spheres. At 
the modest values of S that concern us here, such in- 
stabilities do not occur. It thus remains to study spin- 
odal instabilities of the polydisperse crystal against solid- 
solid demixing. The fact that the solid cloud curve has a 
branch showing solid-solid phase separation already sug- 
gests that such instabilities should be present. Indeed, 
Bartlett found a solid-solid spinodal [2l|. though with a 
thcrmodynamically inconsistent assignment of the excess 
chemical potentials (see Section[nJ . Within the MFE the 
criterion for the spinodal takes its usual form [27j: it is 
the point where the determinant of the curvature matrix 
of the moment free energy, d 2 /„/ (dpidpj) = dpi/dpj, 
first vanishes. The zero eigenvector of the matrix at this 
point gives the instability direction. Using this crite- 
rion, we find the results in Fig. |S](a). The single-phase 
solid is always stable at modest densities or polydisper- 
sities - the spinodal determinant is positive here - but 
can become unstable at larger <f> and 5. With growing 
5, this instability affects a wider and wider range of ip. 
The figure also shows that the spinodal for a triangu- 
lar size distribution is very close to the cloud curve for 
the onset of solid-solid phase separation: past the cloud 
point, a single-phase solid very quickly becomes locally 
unstable. For a Schultz distribution, on the other hand, 
cloud curve and spinodal are well separated as can be 
seen in the inset. This reinforces our above discussion of 
cutoff effects: the latter favour an earlier onset of phase 
separation, cf. Fig. [7| The spinodal condition, on the 
other hand, is known on general grounds (see e.g. pT^) 
to depend only on the moment densities of the parent, 
up to o~ for our excess free energy involving moments 
up to a 3 . Since these parent moments are almost cutoff- 
independent, so is the spinodal curve. (This insensitivity 
of the location of the spinodal is confirmed by the fact 
that the spinodal curves for the triangular and Schultz 
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FIG. 8: Spinodal instability of the polydisperse hard sphere 
crystal against solid-solid demixing, in the volume fraction 
- polydispersity plane {cj>,8). (a) Spinodal (solid) and cloud 
curve (dash-dotted) for triangular (main graph) and Schultz 
(inset) size distributions, (b) The line segments on the spin- 
odal line indicate (for the triangular case) the direction of the 
unstable fluctuations, (c) Comparison of instability direction 
(arrow), path to the "locally optimal" phase (solid line and 
empty circles) , and cloud and shadow solids at the same par- 
ent polydispersity (full circles connected by dotted line). 



cases would be essentially indistinguishable on the scale 
of Fig. 1 (a).) 

We now turn to the nature of the spinodal instability, 
focusing on the case of a triangular size distribution. This 
can be quantified by projecting the instability direction 
at the spinodal into the (</>,5)-plane. The results are indi- 
cated by the line segments on the spinodal line in Fig. [S] 
(b). Bartlett found instability directions which affected 
only the polydispersity 8 while leaving <f> essentially un- 
changed |2l| , which would correspond to vertical lines in 
the plot. By contrast, our analysis shows that the insta- 
bility actually affects both <j> and 5, with relative changes 
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that are of the same order of magnitude. This is con- 
sistent with the properties of coexisting solids discussed 
in Sec. IIV Bl below, which exhibit a strong correlation 
between <f> and S. 

More puzzling is that, in Fig. 0(b), the instability di- 
rections at low S indicate a tendency towards the forma- 
tion of a more polydisperse solid. This appears counter- 
intuitive at first: solid-solid phase separation is driven by 
fractionation and so one expects a preference for smaller 
rather than larger 5. Also, the proximity of the spinodal 
to the cloud curve suggests that the spinodal instability 
direction should be similar to the direction connecting 
the cloud solid and the coexisting shadow. From Fig. 0] 
the instability should therefore point towards larger <j) 
and, again, smaller 5. 

To understand this apparent paradox we consider in 
more detail the "shape" of the MFE f mom at the spin- 
odal, as a function of the moments po, . . . , p$. It is use- 
ful to subtract the tangent plane to /mom at the parent 
phase; the resulting tangent plane distance (tpd) differs 
from / mom only via constant and linear terms in the pi. 
A stable parent is then a local minimum of the tpd, at 
"height" tpd = 0, and any phases coexisting with the 
parent (e.g. the shadow phase(s) for a parent at its cloud 
point) would have the same property. Now, as the spin- 
odal is approached, the curvature of the tpd around the 
parent vanishes in one direction and a "path" towards 
lower, negative, values of the tpd appears; the spinodal 
instability indicates the initial direction of this path. To 
establish where this path leads it makes sense to follow 
it to the nearest "locally optimal" phase, i.e. the nearest 
local minimum of the tpd. If this path is curved in the 
(poj ■ ■ • 5 P3)-space, its initial direction will not necessarily 
capture the properties of the end point, i.e. the locally 
optimal phase. This is the origin of the counter-intuitive 
instability directions that we observe. A specific exam- 
ple is shown in Fig. (c) : the path to the locally optimal 
phase first moves to higher S, consistent with the spin- 
odal instability direction, but the locally optimally phase 
ends up having a smaller polydispersity S than the par- 
ent phase. It also has a larger volume fraction <f>, and the 
change from the unstable parent to the locally optimal 
phase is in a direction comparable to that between cloud 
and shadow, in line with the intuition discussed above. 



B. Phase diagram 

Having clarified the onset of phase separation in poly- 
disperse hard spheres, we next consider the behaviour in- 
side the coexistence region. We have already established 
that, apart from possible cutoff effects, the Schultz and 
triangular parent size distributions give qualitatively sim- 
ilar results, and therefore restrict our attention to the lat- 
ter in the following. Overall features such as the topology 
of the phase diagram should, at the low polydispersities 
5 of interest here, be similar for other size distributions. 

Fig. shows the full phase diagram for the triangu- 
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FIG. 9: Full phase diagram for polydisperse hard spheres with 
a triangular size distribution. In each region the nature of 
the phase(s) coexisting at equilibrium is indicated (F: fluid, 
S: solid). Dashed line: best guess for the phase boundary 
in the region where our numerical data become unreliable. 
From [2^1 . 



lar parent distribution. In each region the nature of 
the phase(s) coexisting at equilibrium is indicated. The 
cloud curves of Fig. 0] (a) reappear as the boundaries 
between single-phase regions and areas of phase coexis- 
tence. Starting from the onset of solid-solid separation 
and increasing density or d, fractionation into multiple 
solids occurs. The overall shape of the phase bound- 
aries in this region is in good qualitative agreement with 
the approximate calculations of [Ifj, see Fig. though 
as discussed below the details of the fractionation be- 
haviour are rather different. We find up to 4 coexisting 
solids. At larger 5 than we can tackle numerically, phase 
splits into 5 or more solids would be expected since each 
individual solid can only tolerate a finite amount of poly- 
dispersity. However, from Fig. [5] such phase splits would 
occur at increasing densities and eventually be limited by 
the physical maximum volume fraction (f> c w 74%. Also, 
at higher 8 more complicated single-phase crystal struc- 
tures, with different lattice sites occupied preferentially 
by (say) smaller and larger spheres, could appear and 
compete with the substitutionally disordered solids we 
consider. 

A feature of the phase diagram in Fig. not predicted 
in previous work is the coexistence of a fluid with multiple 
solids. However, that a three-phase F-S-S region must 
occur was already indicated by the triple point which we 
found earlier on the solid cloud curves. As in the case of 
solid-solid phase splits, coexistences involving more than 
two solids - and a fluid - then appear with increasing S. 

We consider the fractionation behaviour in the mul- 
tiphase regions more systematically in the next section. 
Before doing so, a few qualitative statements arc in or- 
der. In Fig. ED we show a sample plot of the normalised 
diameter distributions n(a) — p(a)/po of four coexisting 
solids. This shows that fractionated solids do not, as one 
might naively assume |l6j |. split the diameter range of 
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FIG. 10: Normalised size distribution of four coexisting solid 
phases obtained from a parent distribution with (<f), 8) — 
(0.63,0.08). From left to right, the solids have volume 
fractions and polydispersities (0.601,0.054), (0.629,0.046), 
(0.646,0.040), (0.663,0.036). From [2(J. 
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FIG. 11: The properties of the daughter solids (circles) aris- 
ing by phase separation from some chosen parents (squares) 
across the phase diagram. Plotted are polydispersity 5 versus 
volume fraction (j>. The arrows show the daughter phases for 
three parents explicitly; as indicated by the dotted lines, not 
all daughter phases are within the range of the plot. Note the 
clustering of all daughter phases near the solid cloud curve. 



the parent evenly among themselves. The polydispersi- 
ties of the coexisting phases are in fact rather different; 
in Fig. 1101 thev range from 8 — 0.036 to 0.054 for a par- 
ent with (5 = 0.08. There is in fact a strong correlation 
between the polydispersity of a fractionated solid and 
its volume fraction: solids with lower volume fraction <p 
tend to have higher polydispersity 8. This conclusion 
is intuitively appealing since higher compression should 
disfavour a polydisperse crystalline packing. 

We have studied the relation between polydispersity 
and volume fraction more quantitatively, by plotting 8 vs 
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FIG. 12: Normalised diameter distributions for F-S-S-S phase 
coexistence obtained from a parent distribution with (<j>, 5) = 
(0.603,0.08). From [H- 



(f> for all the "daughter" solids that arise by phase separa- 
tion from a number of different parents across the phase 
diagram. We find a set of points (Fig. Ill|) that cluster 
very closely around the high-density branch of the solid 
cloud curve, emphasising the tight correlation between 
8 and (j>. Note that some of the points fall above the 
solid cloud curve. This is not a contradiction because 
the latter marks the onset of instability against phase 
separation only for solids with a triangular size distribu- 
tion, whereas the daughter phases plotted here can have 
rather different size distributions (compare Fig. HOfl . 

As part of our qualitative overview of fractionation be- 
haviour, we show next in Fig. 1121 the size distributions 
for a situation where a fluid coexists with three solids. 
The general trend which we observed from the cloud and 
shadow curves, namely for the solid(s) to contain the 
larger particles, is found confirmed here. However, the 
details of the fractionation are again nontrivial: while 
the coexisting fluid is enriched in the smaller particles as 
expected, it also contains "left over" large spheres that 
did not fit comfortably into the solid phases. It thus in 
fact ends up having a larger polydispersity (0.104) than 
the parent (0.08) in this example. 

Finally, an indirect manifestation of fractionation is 
provided by the variation of the osmotic pressure along 
a dilution line. In a monodisperse system, the pressure 
remains constant throughout any phase coexistence re- 
gion because the properties of the coexisting phases do 
not change; only the fractions of system volume vary 
which these phases occupy. In a polydisperse system, on 
the other hand, the composition of the coexisting phases 
varies as the coexistence region is traversed. We illus- 
trate this in Fig. 1131 for a triangular parent size distri- 
bution with 8 — 0.08. It is striking that the variation of 
the pressure with volume fraction is almost smooth, even 
though a number of phase boundaries are crossed. 
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FIG. 13: The osmotic pressure plotted as a function of the 
parent volume fraction along a dilution line, for a triangular 
parent with polydispersity 8 = 0.08. Phase boundaries are 
marked by full circles; line segments are annotated with the 
nature of the phases (fluid/solid) in the different coexistence 
regions. 



V. FRACTIONATION BEHAVIOUR 

We proceed in this section to a systematic study of 
the fractionation behaviour of polydisperse hard spheres, 
having discussed its qualitative features above. To this 
end we extend the classical visual representations in 
terms of cloud and shadow curves and overall phase di- 
agrams to include more detailed information about the 
properties of the coexisting daughter phases. To obtain 
insights into the effects of varying both the parent's vol- 
ume fraction and its polydispersity, 3-D plots will be par- 
ticularly useful here. In a second part we ask whether 
there is an optimal way of making the separation between 
fractionated phases visible, and suggest principal compo- 
nent analysis (PCA) as a method for achieving this. We 
focus throughout on the range of parent polydispersities 
0.04 < S < 0.08, which covers all the various coexistence 
regions in the phase diagram of Fig. [§] Where it is neces- 
sary to distinguish the volume fraction and polydispersity 
of the parent from those of the daughter phases, we will 
add the superscript (0), writing </>(°) and S^°\ 

We start with a 2-D plot showing the volume fraction 
of the coexisting phases versus the volume fraction of the 
parent phase, Fig. FLU (a), for two parent polydispersities 
8. For a narrow parent size distribution (8 — 0.04, inset), 
we see that the behaviour in the F-S coexistence region 
is similar to what would be expected for a monodisperse 
system, with the volume fractions of the daughter phases 
remaining essentially constant. Only at large <f>( ' does 
the polydisperse nature of the system become fully appar- 
ent, through the occurrence of S-S phase separation. For 
a parent with 8 = 0.08, on the other hand, the properties 
of the daughter phases vary strongly with 0'°' . In the 
S+S+S and S+S+S+S regions in particular, the volume 
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FIG. 14: (a) Volume fraction (j> of the various coexisting 
daughter phases versus the volume fraction (j>^ of the parent 
phase, for 8 = 0.08 (main graph) and 5 — 0.04 (inset), (b) 3-D 
plot showing the dependence of the (j> on 4>^ and the parent's 
polydispersity 8. Different phases are represented by different 
grey levels. Note that the 5-axis is plotted upside down for 
better visibility. The top and bottom slices correspond to the 
2-D plots in (a). 



fractions of the daughter solids increase systematically 
with </>(°) : fractionation from a denser parent here pro- 
duces denser daughter phases, rather than varying pro- 
portions of daughters with fixed densities. 

To demonstrate more explicitly the change in be- 
haviour as the parent polydispersity 8 increases, we show 
in Fig. [21 (b) a 3-D plot of the daughter volume fractions 
4> versus and 5. The orientation of the axes has been 
chosen such that horizontal cuts through the plot repre- 
sent fixed 8, with the top and bottom planes correspond- 
ing to the data shown in the 2-D plots of Fig. ^] (a). 
A benefit of the 3-D representation is that each daugh- 
ter phase now corresponds simply to a separate surface. 
Each surface ends at the phase boundary where the rel- 
evant phase disappears from the phase split. The disap- 
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FIG. 15: Mean diameter of coexisting phases plotted against 
parent volume fraction (fi^, for parent polydispersity 8 — 
0.08. The dashed lines delineate the various phase coexistence 
regions. 



pearance or appearance of any phase then causes kinks 
in the other surfaces. As expected, only the fluid sur- 
face extends to the lowest 4>^ ■ As 4>^ is increased, the 
"conventional" solid which also exists in the monodis- 
perse limit makes its first appearance. A further three 
fractionated solids then eventually appear one after the 
other. These are polydispersity-induced, i.e. have no ana- 
logue in the monodisperse system, and the surfaces rep- 
resenting them do not extend to 8 — > 0. 

Having clarified the variation of the volume fractions of 
the daughter phases across the phase diagram, we show 
their mean diameters in Fig. El plotted against parent 
volume fraction at fixed (parent) polydispersity S = 0.08. 
One observes clearly the general trend for the solid phases 
to contain larger particles than the fluid. An exception 
to this occurs in the F+S+S+S coexistence region, where 
the fluid has a slightly larger mean diameter than one of 
the solids. The explanation for this can be found in our 
earlier discussion of Fig. El in addition to the small- 
est spheres, the fluid can also contain some of the larger 
spheres that are not accommodated in any of the solids, 
and this pushes up its mean diameter. The second quali- 
tative trend demonstrated by Fig. El is that the coexist- 
ing solids tend to split the range of particle diameters in 
the parent distribution amongst themselves, with almost 
equidistant mean diameters. As the parent volume frac- 
tion increases, the strength of this fractionation effect is 
seen to grow, and the mean diameters become increas- 
ingly separated from each other. 

Finally, we examine the relationship between the poly- 
dispersities 6 of the different daughter phases and the 
volume fraction </> ( - - ) and polydispersity 5^ of the par- 
ent. Fig. El (a) shows 2-D plots of the daughter poly- 
dispersities versus 4>^°\ at 5^ — 0.04 and 0.06. As ex- 
pected, for the more polydisperse parent there are signifi- 
cant variations of the daughter polydispersities across the 
coexistence regions. Where multiple solids coexist, their 




FIG. 16: (a) Plot of the polydispersities of coexisting phases 
along two dilution lines (thin horizontal lines), i.e. as a func- 
tion of the parent volume fraction for fixed parent polydis- 
persity <5<°) = 0.04 and 0.06. The dashed lines indicate the 
phase boundaries in the (</> , <r°')-plane; phases appear or 
disappear at the points where the horizontal line correspond- 
ing to the fixed parent polydispersity intersects these phase 
boundaries (full circles), (b) Corresponding 3-D plot, show- 
ing the daughter polydispersities 5 against the parent volume 
fraction cjy- ' and parent polydispersity 8^°\ 



polydispersities decrease with increasing parent volume 
fraction. This is consistent with the general trend that 
denser solids tend to be less polydisperse. In the 3-D plot 
of Fig. mSl fright), this same trend also causes the surfaces 
corresponding to the various solids to have rather similar 
shapes in the region of solid-solid coexistence. For the 
fluid, on the other hand, the graph demonstrates that it 
always has a larger polydispersity than the parent, aris- 
ing from the presence of large particles "left over" from 
the solid phases (see Fig. I12|) . 
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A. Principal components analysis 

The above plots of aspects of fractionation behaviour 
lead naturally to the question of whether there is a 
"maximally fractionating" property, i.e. one which most 
strongly reveals the differences between the various co- 
existing phases across the phase diagram. We focus on 
properties which are generalised moments of the density 
distribution, of the form r = Jda f(a)p(a) with some 
weight function /(cr). While not all properties can be 
expressed in this way - the polydispersity 5, for example, 
involves squares and ratios of such moments - this is still 
a fairly large class of measurable properties; e.g. setting 
/(cr) = 1 would give us the number density, /(cr) = cr 3 
the volume fraction, /(cr) = a the mean diameter times 
the number density etc. 

Suppose now that we have a number of measurements 
of p(cr), specifically the daughter density distributions 
that arise within some region of the phase diagram. We 
can think of the p(a) as points in a high-dimensional (in 
fact infinite-dimensional) space, and of our desired mo- 
ment r as a projection along the direction defined by 
/(cr) [2?J. A good choice for a maximally fractionating 
property would then be to maximise the variance of our 
moment among the various measured p(a). This can be 
done by Principal Component Analysis (PCAla method 
designed to select directions of large variance [4!| . Math- 
ematically, the requirement of maximum variance can be 
written as 

max /(cr) J da da' f(a)A(a,a')f(a') (18) 

subject to Jdaf 2 (a) — 1; here A(a,a') is the (infinite- 
dimensional) covariance matrix of our measurements. We 
define this as A(a,a') = ([p(a)-p^{a)][p(a')~p^{a')}). 
The average here is over all our measurements of p{a), 
and we subtract off for each p{a) the corresponding par- 
ent distribution p^ (a) . This effectively removes the av- 
erage of the various measured p(a) because, from particle 
conservation JSJl, the parent is a weighted average over 
the various daughter phases. An alternative definition of 
A would be to remove the actual measurement average, 
A(a : a>) = {[p(a) - (p(a))} [p(a>) - (p(a'))}). In our nu- 
merical experiments described below, this lead to almost 
indistinguishable results. 

The maximisation problem (|18fl is in principle over 
an infinite-dimensional function space. To arrive at a 
more practical task, we restrict the search to a sub- 
space by requiring the weight function to be of the form 
f( a ) — Si=o ai<Tl - This corresponds to searching for a 
maximally fractionating property among those express- 
ible as linear combinations of the moments po, . . . , p3, i.e. 
r = X)j=o a iPi With this simplification, the problem 118|) 
reduces to 

max a a T Ca subject to a 1 Da =1 (19) 
Here a denotes the vector with elements olq, . . . , 0:3 and 
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FIG. 17: Principal component weight function, f(a), obtained 
from data along dilution lines for three different values of 
parent polydispersity 0.04, 0.06 and 0.08. The range of a- 
values where f(u) 7^ is in each case that over which the 
parent density distribution is nonzero, i.e. the range of particle 
sizes actually occurring in the system. 



the 4x4 matrix C is defined as 

Cy = Jda da' a*A(a 7 a'){a'y = (( Pi - p^)( Pj - pf)) 

which is just the covariance matrix of the moments, with 
the parent moments again subtracted off. The matrix 
D, on the other hand, is given by Dy — J da a' l+ K The 
cr-integration range has to be bounded to make this well- 
defined. In our case of a triangular parent distribution 
the obvious choice, adopted here, is to make this range 
equal to the range of particle sizes occurring in the par- 
ent. 

Imposing the constraint in l|19f) via a Lagrange multi- 
plier shows that solution vectors a must obey Ca — XDa, 
or equivalently D~ 1 l 2 CD~ 1 / 2 {D 1 / 2 a) = \D 1 / 2 a. The 
solutions can thus be obtained by an eigenvalue decom- 
position of the matrix D~ 1 / 2 CD~ 1 / 2 , with A the eigen- 
value and D 1 / 2 a the corresponding eigenvector. (Nu- 
merically, it is more convenient to solve the equivalent 
problem of finding the eigenvalues and right eigenvectors 
of the matrix D~ 1 C.) The eigenvectors are termed prin- 
cipal components, and the A's give the variance captured 
by each principal component. The most important prin- 
cipal component, and the one of interest to us, is then 
the one with the largest A. 

We have implemented this PC A search for maximally 
fractionating properties by considering as our measured 
p(a) the daughter phases as they occur along a dilution 
line. We do this separately for triangular parent distribu- 
tions of polydispersity 0.04, 0.06 and 0.08, respectively, 
because different ranges of particle size a are relevant 
in the three cases. The resulting weight functions /(cr) 
are plotted in Fig. El O ne sees that in all cases, f(a) 
is to a good approximation a combination of the odd 
weight functions a and cr 3 , with the coefficients such that 
/(cr) crosses zero near the edge of the cr-range. Loosely 
speaking, the function /(cr) can be interpreted as an ap- 
proximation to sgn(cr — 1) within the space spanned by 
cr°, . . . , cr 3 , i.e. by a third-order polynomial in a. It thus 
effectively measures the difference in number density be- 
tween particles above and below the mean parental di- 
ameter. This is an intuitively appealing measure of frac- 
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FIG. 18: Maximally fractionating moment (as selected by 
PCA) for coexisting daughter phases, plotted against parent 
volume fraction (j>^ at parent polydispersity S = 0.04, 0.06 
and 0.08 (from left to right). 



tionation behaviour. 

Finally, Fig. ^3] shows the properties of the daughter 
phases as measured by the maximally fractionating ob- 
servable selected by PCA. The overall features of the plot 
on the right, for parent polydispersity 6 — 0.08, are not 
dissimilar to the mean diameter representation in Fig. 1151 
so that the benefit of PCA in this problem is relatively 
modest. Some interesting features are accentuated by 
PCA, however; e.g. the crossover between the fluid and 
solid lines is more pronounced in Fig. 1181 demonstrating 
clearly how the fluid size distribution acquires a signifi- 
cant fraction of the larger particles. We expect that the 
benefits the PCA method of selecting maximally frac- 
tionating properties should become more pronounced in 
systems with several polydisperse attributes ct, e.g. parti- 
cle size and charge. Suitable properties for revealing frac- 
tionation behaviour could then depend on combinations 
of these attributes, which can be systematically found 
using PCA. 



VI. COMPARISON WITH PERTURB ATIVE 
THEORIES AND MONTE CARLO 
SIMULATIONS 

In this section we validate our theoretical predictions in 
two ways. First, we compare to perturbative theories for 
near-monodisperse parents, which predict qualitatively 
how the properties of coexisting phases should vary as 
the parent is made increasingly polydisperse. Second, 
we compare quantitatively to Monte Carlo simulations 
of polydisperse hard spheres with an imposed chemical 
potential distribution. 



A. Near-monodisperse systems 

For systems that are nearly monodisperse, one can 
make general statements about the fractionation be- 
haviour by considering deviations of particle diameters 
from the mean as small and performing a perturbation 
expansion [50ll5l|. This approach presupposes, of course, 
that the phase separation of interest occurs already in the 
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FIG. 19: Log-log plots of the difference in mean diameters 
(left) and in polydispersity (right) for coexisting fluid-solid 
phases at parent density p = 0.52, plotted against parent 
polydispersity 8. The solid lines show the theoretically ex- 
pected power laws, Eqs. I20H and 12 U . with the proportional- 
ity constants fitted by eye to the numerical data. The insets 
show the ratios Aa/S 2 , AS/S 3 , which from the theory are 
expected to approach constants for S — > 0. The data are 
consistent with this (note the narrow ranges displayed on the 
y-axes). 



monodisperse reference system. In our hard sphere case 
it is therefore applicable only to fluid-solid coexistence. 
Phase separation involving several solid phases is induced 
by polydispersity itself and cannot be treated perturba- 
tively. 

A strong prediction of the perturbative approach is 
that the difference in mean particle diameters of two co- 
existing phases, Act = — is universally related 
to the parental polydispersity 5 via 



Act 



(20) 



An increase in the width of the parent size distribution 
thus contributes only at second order to the mean di- 
ameter difference. The proportionality coefficient in this 
relation is non-universal and depends on the properties 
of the phase separation in the corresponding monodis- 
perse reference system. Equation (|20[) is the leading term 
in a perturbation expansion in 6 at fixed density of the 
parent. The volume fraction (j) then varies with 5; for 
symmetric size distributions it increases. If instead <fi is 
held constant as S is varied, the perturbation expansion 
is modified, though the leading term (|20|1 remains unaf- 
fected. 

A relation analogous to equation (I20fl applies to the dif- 
ference in polydispersities A<5 = S^ 1 ' — S^ of the daugh- 
ter phases. For symmetric parent size distributions, one 
finds 



AS 



(21) 
only 



so that an increase in the parent polydispersity S 
affects AS at third order. 

To verify the above perturbative predictions, we de- 
termined the mean diameters and polydispersities of the 
daughter phases for fluid-solid separation at parent den- 
sity p = 0.52 and parent polydispersity S ranging from 
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FIG. 20: Plot of the polydispersity of coexisting solids versus 
the polydispersity 8 of the parent phase. The parent number 
density is kept constant at p = 0.62 while the latter is varied. 
The inset shows the log-log plot of the difference of the poly- 
dispersities AS in the S+S and S+S+S region versus S and 
the slope of the theoretically expected power law 121H (dashed 
line). This demonstrates that the perturbative prediction for 
the scaling of AS with S does not apply here. 



0.02 to 0.086. Fig. m shows log-log plots of Act and AS 
against S, confirming the predicted power laws. The in- 
sets show the ratios Aa/S 2 and AS/S 3 . Our data are 
again consistent with the theoretical expectation that 
these ratios should approach constants in the limit of 
small 5. Overall, the scaling predictions of perturbative 
theories for near-monodisperse systems are fully obeyed 
by our data. 

Finally, to emphasise the point that perturbative ap- 
proaches do not apply to phase separations that are 
caused by polydispersity, we show in Fig. [201 the evo- 
lution of the polydispersity of coexisting solids as the 
parent polydispersity is increased, at fixed parent den- 
sity p = 0.62. As the inset shows, there is now no simple 
relationship akin to (|21|l which relates the difference in 
the polydispersities of the daughter solids to the S of the 
parent. In fact, the perturbative limit of small 5 is not 
even defined here, since the solid-solid phase separation 
only occurs above a nonzero (density-dependent) thresh- 
old value of 5. 



B. Comparison with Monte Carlo simulations 

As discussed in the preceding sections, our theoreti- 
cal results for fluid-solid coexistence in polydisperse hard 
sphere systems are in qualitative agreement with numer- 
ical simulations 0, 0, 1 1 1 1] , in particular concerning the 
coexistence of rather polydisperse fluids with solids that 
have a much narrower size distributions. There is, how- 
ever, an important difference: our calculations apply to 
the experimentally realistic case where an overall par- 
ent density distribution is fixed. The simulations, on 



the other hand, are performed at imposed chemical po- 
tential differences, with the actual size distributions in 
the coexisting phases varying strongly across the phase 
diagram. In order to obtain a quantitative comparison 
between theory and simulations, we calculate explicitly 
in this section the theoretical predictions for the - some- 
what unrealistic - scenario addressed in the simulations. 
We will find good quantitative agreement, thus validat- 
ing our approach and, in particular, our choice of free 
energy expressions for the fluid and solid phases. 

The simulations of 0, 0, |n| were carried out in an iso- 
baric semi-grandcanonical ensemble, which corresponds 
to fixed particle number N, pressure P and chemical po- 
tential differences fx(u) — /z(ct{,). Here a b is the diame- 
ter of a reference particle. The advantage of the semi- 
grandcanonical ensemble is that it allows many different 
realizations of the particle size distribution to be sam- 
pled, thus minimising finite-size effects. The fixed par- 
ticle number N, on the other hand, avoids simulation 
moves where particles need to be inserted into dense flu- 
ids or solids. 

Bolhuis and Kofke 0, Q considered specifically a 
quadratic form for the chemical potential differences, 



2v 



(22) 



The activity exp[/i(er)] thus has a Gaussian shape of vari- 
ance v. For small v, one expects the activity distribu- 
tion to set the size distributions in the coexisting phases, 
which should therefore have polydispersity S — v 1 ' 2 ] 
v — > recovers the monodisperse case. The reference 
diameter a b = 1 is held fixed as v is increased from zero. 
The pressure P is then adapted by Gibbs-Duhem integra- 
tion [52| to follow the line of fluid-solid phase coexistence 
in the (u, P)-plane. 

In order to reproduce the situation considered in the 
simulations using our theoretical approach, we will study 
a system with prior R(a) — exp[— (<r — I) 2 /(2v)\. The 
moment free energy then gives the free energy of phases 
with density distributions of the form (cf. Q12[l) 



p(a) = cxp 



(^-1) S 
2v 



J>ct< 



i=0 



From © , the corresponding chemical potentials have the 
form 



jix(a-) = lnp(a) ■ 



»=0 



(ct-1) 



2v 



Now the Xi + /zf x are just the moment chemical poten- 
tials pi = df mom /dpi. So if we apply the moment free 
energy but treat the moment densities pi, p2, P3 as non- 
conserved, the associated pii are forced to vanish auto- 
matically at equilibrium. The sum over i in (|23|l then re- 
duces to a constant, po — Xq +Pq x , and we have precisely 
the chemical potential differences l|22|) used in the sim- 
ulations. In summary, applying the MFE method with 
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a Gaussian prior and po the only conserved moment, we 
increase po from zero until coexistence with a solid phase 
is first found. This is then the desired fluid-solid coex- 
istence for quadratic chemical potential differences, and 
we can determine in particular the pressure P at coexis- 
tence. Repeating this process for a range of values of v 
gives the coexistence curve in the (y, P)-plane. 

Our actual implementation of this approach has one 
minor difference. In the simulations it is observed that 
the mean diameters in the coexisting phases decrease sig- 
nificantly as v is increased, eventually becoming much 
smaller than ov For our numerical work, however, it 
is desirable to keep the size distributions within a fixed 
range, e.g. in order to ensure that our chemical poten- 
tials for the solid phase remain reliable. To achieve this, 
we treat not just po but also pi as conserved. Keeping 
Pi I Po — 1 as Po is varied then ensures that the fluid 
phase always has unit mean diameter, and the particle 
sizes in the coexisting solid are expected to be compara- 
ble. This ensures that we can use a fixed er-range for all 
calculations, for which we choose a £ [0.7, 1.3]. 

With po and p\ both conserved, the chemical poten- 
tials iffi^ become 

= ~ o +P0+PKJ (24) 

(a — 1 — vi±\ ) 2 1 o , 
= 2y — +-^MI+Mi+Mo (25) 

which is again of the form 1)221) but now with a vary- 
ing reference diameter 05 = 1 + vpi- The corresponding 
scaled quantities that are to be compared to v and P from 
the simulations are then vja\ and Pa\ [5{|. Note finally 
that our numerical implementation again uses centred 
moments, with weight functions [{a — l)/<5o] 4 rather than 
cr% but this causes no conceptual differences. In particu- 
lar, keeping the standard moments po and p\ conserved is 
equivalent to conservation of the centred moments with 
i = and i = 1, because of the linear relations between 
the two sets of moments. 

Fig. (a) shows our results for the coexistence curve 
in the (yja\ , Pu^-plane. As v increases (starting from 
the bottom left corner), both Per 2 and v ja\ initially in- 
crease. However, eventually vja 2 reaches a maximum 
value ^max = V I&1 = 0.0056. At this point, the slope 
d(Pa%) J ' d(y / 'a 2 ) becomes infinite. On further increas- 
ing v, the coexistence curve then bends back, with vja 2 
decreasing towards zero while the pressure diverges. Bol- 
huis and Koike 8] argued that this divergence arises be- 
cause the pressure is measured on the scale of the mean 
at, of the activity distribution, while the typical particle 
diameters in the coexisting phases become much smaller 
than <Jb, by a factor scaling as v /a 2 . The rescaled pres- 
sure Paf (i//ct 2 ) 3 — Pv 3 ja\ should therefore approach a 
constant value in the limit v/<7 2 — > 0. The simulations 
were consistent with this expectation, and our theoretical 
results in Fig. ED(b) are m full agreement. By extrapo- 
lation, we estimate the limiting or 'terminal' value of the 
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FIG. 21: (a) Solid-fluid coexistence pressure Po\ as a function 
of the imposed width v/a^ of the activity distribution. Both 
are scaled appropriately with at to account for the fact that 
Ob varies in our calculation but is held constant in the simula- 
tions, (b) The pressure is rescaled to Po\(y / 'of) 3 = Pf 3 /(jf 
to show the limiting behaviour for v/cr% — > 0. 



rescaled pressure, i.e. the point where the rescaled coexis- 
tence curve intersects the vertical axis, asPt = 7.9xl0~ 5 . 

The scaling mentioned above implies that, as v be- 
comes large, the mean particle diameter in the coexist- 
ing phases will be of order 06(^/af) = v/cxf,, rather than 
(7b- I n our scheme, where the mean diameter in the fluid 
is fixed at unity, <7b should thus become linear in v. As 
shown in Fig. [53 (left), this is indeed what we find. A 
plot of the numerical derivative of this dependence, in 
the inset of Fig. [221 (left), also reveals that at small re- 
values - below those where v ja\ reaches its maximum - 
the behaviour is no longer exactly linear. This is to be 
expected considering that (7b = 1 + vp\ depends on both 
v and Hi. 

The plots in the middle and on the right of Fig.l22lshow 
particle size distributions in the coexisting fluid and solid 
phases at two different values of v. For small v = 0.02 2 , 
the distributions are essentially identical and have width 
S w v 1 / 2 = 0.02 as expected; they are also close to the 
activity distribution, which has its peak at <7b = 1-02 
for this v. For larger v = 0.08 2 , on the other hand, 
there is significant fractionation between the fluid and 
the solid. One can now also clearly see how the mean 
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FIG. 22: Left: Plot of at, as a function on v. The inset shows 
the derivative do^jdv. As expected, at grows linearly as v 
becomes large. The star indicates the value of v at which 
v /of reaches its maximum and the slope of the pressure plots 
in Fig. l2"Tl becomes infinite. Middle and right: Normalised size 
distributions n(a) of the coexisting fluid (dashed line) and 
solid (solid line) phases. The dotted curve gives the shape of 
the activity distribution exp(/x(cr)). Middle: For v — 0.02 2 , 
fluid and solid have essentially identical size distributions of 
polydispersity 8 ~ v 1 ^ 2 = 0.02; the corresponding value of 
fj b is 1.02. Right: For v — 0.08 2 , the size distributions are 
significantly different from each other and from the activity 
distribution, which is now centred around db = 1.19. 
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FIG. 23: Phase diagram for fluid-solid coexistence with im- 
posed quadratic chemical potential differences. Plotted is the 
polydispersity 8 versus the volume fraction <j> of the coexist- 
ing fluid and solid phases; v increases from bottom left to 
top right along the curves. The circles indicate the terminal 
points reached by extrapolating to v — * oo. The dotted lines 
sketch the approach to the known monodisperse limit v — > 0. 



particle diameters - which are exactly unity in the fluid, 
by construction, and around 1.02 in the solid - become 
smaller than the mean of the activity distribution, which 
is (Tfe = 1.19 for this value of v. 

To summarise the properties of the coexisting fluid 
and solid phases, we plot them in a volume fraction- 
polydispersity phase diagram, shown in Fig. 1231 As dis- 
cussed in detail in Q , a feature which is at first surprising 
is that the curves terminate, with the properties of the 
coexisting phases approaching finite limits for v — > oo. 
One has to bear in mind, however, that the occurrence 
of such terminal points is directly linked to the shape of 



Quantity 


Bolhuis and Kofke [8] 


Present work 


"max 

Pt 

(</Vs, 5t, s ) 


0.0056 
7.9 xlO" 5 
(0.545, 0.12) 
(0.575, 0.057) 


0.0056 
7.9 xlO" 5 
(0.548, 0.113) 
(0.592, 0.057) 



TABLE II: Comparison between some characteristic quan- 
tities of fluid-solid coexistence at imposed chemical poten- 
tial distribution, as determined in simulations |£| and in the 
present theoretical study. Here f ma x is the maximum value of 
via 2 , P t the terminal value of the rescaled pressure Pv 3 ja\ 
in the limit v/a 2 — > 0, 4>t,t/s the terminal volume fraction for 
the fluid/solid phases and #t,f/s the corresponding terminal 
polydispersity. 



the imposed chemical potential distribution and so not 
physically very meaningful. Indeed, other shapes can and 
do give fluids and solids with larger </> and/or S 

We obtain the location of the terminal points by plot- 
ting our numerical predictions against 1/v and extrapo- 
lating to \ jv = 0. The resulting values are compared in 
Tablelnlwith those obtained in the simulations of @. We 
find excellent quantitative agreement for ^ max , defined as 
the maximum value of v/cr^, and Pt, the terminal value 
of the rescaled pressure Pv z ja\. Similar comments ap- 
ply to the volume fractions and polydispersities at the 
terminal points of the fluid and solid coexistence curves 
in Fig. [23 Only the terminal volume fraction of the solid 
is over-estimated somewhat, but even here the deviation 
is less than 3%. 

In conclusion, our theoretical predictions for fluid-solid 
coexistence at imposed chemical potential differences are 
not just in qualitative but in fact quantitative agreement 
with the outcomes of Monte Carlo simulations. This pro- 
vides strong validation for our approach. It demonstrates 
in particular that our chosen model free energies for the 
hard sphere fluid and solid are accurate, at least in the 
range of relatively small polydispersities studied here. 

VII. CONCLUSION AND OUTLOOK 

We have studied the equilibrium behaviour of size- 
polydisperse hard spheres, starting from accurate free 
energy expressions for the hard sphere solid and fluid. 
Cloud and shadow curves, which locate the onset of phase 
coexistence, were found exactly by using the moment free 
energy (MFE) method. We were also able to calculate 
the full phase diagram, however, by using the MFE re- 
sults as starting points for a solution of the full phase 
equilibrium equations. 

In contrast to earlier simplified theoretical treatments, 
we found no point of equal concentration between fluid 
and solid. Rather, the fluid cloud curve continues to 
larger polydispersities while the coexisting solid shadow 
always has a polydispersity S below a "terminal" value of 
around S t w 0.06. In this sense the concept of terminal 
polydispersity only applies to the solid phase, while any 
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experimentally observed terminal polydispersity from the 
fluid side must be attributed to non-equilibrium effects 
such as an intervening kinetic glass transition , large 
nucleation barriers 23] or the unusual growth kinetics of 
polydisperse crystals |2J|. 

Concomitant with the absence of the point of equal 
concentration, we also found no re-entrant melting. In- 
stead, a sufficiently compressed polydisperse solid frac- 
tionates into two or more solid phases; our results in this 
region of the phase diagram are consistent with previous 
approximate calculations. In addition, we found that co- 
existence of several solids with a fluid phase is also pos- 
sible. That such phase splits must exist is clear from the 
fact that the solid cloud curve has two branches, describ- 
ing onset of fluid-solid and solid-solid phase separation 
at low and high densities respectively; a fluid-solid-solid 
coexistence region begins where these meet. 

We then analysed the fractionation behaviour in detail. 
As a general rule, the fluid phases contain the smaller 
particles in the system, while the larger ones are found 
predominantly in the solid phases. The solid phases have 
smaller polydispersities 6 than the parent phase; this is 
as expected since narrower particle size distributions are 
more easily accommodated on a regular lattice. Consis- 
tent with this physical intuition, we also found that there 
is a strong correlation between the polydispersity S and 
the volume fraction <f> of coexisting solids, with the denser 
phases (larger cf>) having smaller 5. For the fluid phases, 
on the other hand, we found larger polydispersities than 
in the parent. This is because the fluid contains, together 
with a relatively narrow distribution of smaller particles, 
also residual larger particles that were not incorporated 
into any of the solid phases. 

Three-dimensional fractionation plots transparently 
showed the continuity of the properties of the various 
phases across the phase diagram, with each correspond- 
ing to a distinct surface. The individual phases change 
significantly as coexistence regions are traversed; this 
is in contrast to monodisperse systems, where only the 
amounts of coexisting phases vary. Correspondingly, the 
pressure in the polydisperse case was seen to vary al- 
most smoothly on traversing several coexistence regions, 
whereas it would be constant within each for a monodis- 
perse system. We finally proposed a method for con- 
structing maximally fractionating observables, i.e. mea- 
surable properties which reveal most clearly the differ- 
ences between the various coexisting phases. This was 
based on Principal Components Analysis in the space of 
the relevant density distributions. The benefits of this 
method were modest in our case, but it could be of sig- 
nificant interest for analysing systematically the phase 



behaviour of systems with more than one polydisperse 
attribute, e.g. particle size and charge. 

In the final section we compared our predictions to per- 
turbative theories for near-monodisperse systems, finding 
full agreement. We also performed a detailed compari- 
son with Monte Carlo simulation carried out at imposed 
chemical potential distribution, where particle size dis- 
tributions vary across the phase diagram. The excellent 
agreement obtained provided strong validation of our ap- 
proach and in particular of our choice of model free en- 
ergies for polydisperse hard sphere fluids and solids. 

There are a number of possibilities for extending and 
complementing the present work. Our study was limited 
to systems with relatively narrow size distributions, with 
polydispersities 8 up to « 0.14. At higher i5, fluid-fluid 
demixing would eventually be expected to occur [43, |48| . 
So far only the spinodals for this have been calculated, 
however, and it would be interesting to understand the 
topology of the full phase diagram in this large-^ region. 
One might, for example, expect to find coexistence of 
multiple fluids, but the conditions required for this are 
at present unclear. 

Quantitative studies of the phase behaviour of hard 
spheres at large S would require accurate model free en- 
ergies for wide particle size distributions. For the fluid, 
the BMCSL approximation may continue to be sufficient, 
although a recent comparison with simulations has re- 
vealed some shortcomings |53| . Much more pressing is 
the need for an accurate free energy for strongly polydis- 
perse hard sphere solids. This would allow one to investi- 
gate, for example, whether the dominance of the largest 
particles at the onset of solid-solid coexistence which we 
found for Schultz size distributions is a genuine physical 
effect. A quantitative verification of the prediction that 
polydisperse hard spheres with sufficiently fat-tailed size 
distributions split off multiple fractionated solids even at 
low density [54( would also be of interest. A significant 
challenge in the construction of approximate free energies 
for hard sphere solids is that the simplifying assumption 
of a substitutionally disordered structure - which was 
implicit in our study - may break down at large polydis- 
persities. Competing substitutionally ordered structures 
would then also have to be considered. 

Finally, it will be exciting to generalise our approach to 
more complex colloidal systems, by for example including 
attractive interactions or extending the scope to polydis- 
perse colloid-polymer mixtures. Work on these scenarios 
is currently under way. 

We are grateful to Paul Bartlett for making his numer- 
ical code for the solid free energy available to us, and to 
EPSRC for financial support (GR/R52121/01). 
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